ARTICLE  IN  PRESS 


Combustion  and  Flame  xxx  (2010)  xxx-xxx 


ELSEVIER 


Contents  lists  available  at  ScienceDirect 

Combustion  and  Flame 

journal  homepage:  www.elsevier.com/locate/combustflame 


Combustion 
and  Flame 


Simulations  of  flame  acceleration  and  deflagration-to-detonation  transitions  in 
methane-air  systems 

D.A.  Kessler*,  V.N.  Gamezo,  E.S.  Oran 

Laboratory  for  Computational  Physics  and  Fluid  Dynamics,  Naval  Research  Laboratory,  Washington,  DC,  United  States 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  2  February  2010 

Received  in  revised  form  17  March  2010 

Accepted  16  April  2010 

Available  online  xxxx 


Keywords: 

Deflagration-to-detonation  transition 
Methane-air  explosions 
Flame  acceleration 
Obstructed  channels 


Flame  acceleration  and  deflagration-to-detonation  transitions  (DDT)  in  large  obstructed  channels  filled 
with  a  stoichiometric  methane-air  mixture  are  simulated  using  a  single-step  reaction  mechanism.  The 
reaction  parameters  are  calibrated  using  known  velocities  and  length  scales  of  laminar  flames  and  deto¬ 
nations.  Calculations  of  the  flame  dynamics  and  DDT  in  channels  with  obstacles  are  compared  to  previ¬ 
ously  reported  experimental  data.  The  results  obtained  using  the  simple  reaction  model  qualitatively,  and 
in  many  cases,  quantitatively  match  the  experiments  and  are  found  to  be  largely  insensitive  to  small  vari¬ 
ations  in  model  parameters. 
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1.  Introduction 

The  atmospheres  of  confined  regions  in  underground  facilities, 
such  as  sealed-off  tunnels  and  chambers  in  mining  operations, 
can  develop  into  potentially  explosive  mixtures  of  natural  gas 
and  air.  Explosions  in  these  regions  are  a  significant  concern  be¬ 
cause  of  the  extent  to  which  they  may  harm  personnel,  equipment, 
and  the  production  process.  This  paper  describes  the  first  steps  we 
have  taken  in  developing  a  multidimensional  numerical  model  to 
study  explosions  in  large-scale  systems  containing  mixtures  of  nat¬ 
ural  gas  and  air.  This  final  model  must  have  the  ability  to  compute 
the  different  stages  of  the  evolution  of  a  chemically  reactive  flow: 
ignition  by  a  small  spark,  rapid  flame  acceleration,  development  of 
shocks,  shock-flame  interactions,  and  detonation  initiation  in  com¬ 
plex  geometries  for  natural  gas  mixtures  with  spatially  and  tempo¬ 
rally-varying  stoichiometries.  Here  we  describe  the  first  few  steps 
toward  achieving  this  objective. 

Numerical  models  that  can  describe  the  behavior  of  shocks  and 
detonations  vary  widely  in  their  complexity,  but  for  many  practical 
situations,  an  extensive  description  of  the  details  of  the  chemical 
pathways  is  unnecessary.  Instead,  it  is  more  important  to  have 
an  accurate  model  of  the  fluid  dynamics  coupled  to  a  model  for 
the  chemical-energy  release  that  puts  the  released  energy  in  the 
“right”  place  in  the  flow  at  the  “right”  time.  For  example,  ignition 
behind  a  shock  forming  a  detonation  wave  can  be  quantitatively 
predicted  if  an  acceptable  representation  of  the  chemical  induction 
time  is  known  as  a  function  of  the  state  variables,  temperature,  and 
pressure.  This  observation  led  to  the  development  of  single-step 
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reaction  models,  such  as  the  induction-parameter  model  [1],  and, 
later,  to  reduced  two-  and  three-step  chain-branching  reaction 
models  [2].  While  these  models  have  been  quite  successful  for 
computing  steady-state  and  certain  transient  properties  of  detona¬ 
tions,  they  generally  cannot  be  used  for  calculating  properties  of 
flames  for  which  diffusion  and  thermal  conduction  are  important. 
They  are  not  appropriate  for  combustion  wave  transitions,  such  as 
the  transition  from  a  laminar  to  turbulent  flame  or  a  turbulent 
flame  to  a  detonation. 

Thus,  numerical  simulations  involving  flames  require  some 
treatment  of  diffusion  processes  in  addition  to  chemical  reaction 
and  energy  release.  This  can  be  done  with  high  accuracy  using  de¬ 
tailed  chemical  reaction  mechanisms,  if  they  are  known,  although 
the  computational  price  can  be  so  high  that  it  severely  limits  the 
extent  of  a  calculation.  For  this  reason,  it  becomes  much  too  expen¬ 
sive  to  use  a  detailed  reaction  mechanism  to  compute  flames  in 
large  physical  systems.  An  enormous  amount  of  work  has  been 
done  to  find  reduced  chemical  models  for  hydrogen  and  hydrocar¬ 
bons  that  work  well  when  coupled  to  fluid  codes.  Usually,  this 
amounts  to  finding  the  minimal  set  of  reactions  that  are  critical 
for  describing  the  flame  structure  for  a  particular  fuel  mixture 
[3-6].  A  different  approach  is  to  assume  that  the  complex  set  of 
reactions  can  be  modeled  by  a  generic  set  of  global  reactions  [7- 
9].  For  example,  the  widely  used  model  developed  by  Westbrook 
and  Dryer  [7]  is  based  on  a  single-step  Arrhenius  reaction, 

Q  =  /\exp(-£a/RT)[Fuel]“[Oxidizeif  (1) 

where  [m]  represents  the  concentration  of  species  m.  The  model 
parameters,  A,  Ea,  a,  and  b,  are  calibrated  for  a  particular  fuel  based 
on  measured  lean  and  rich  explosion  limits  and  the  laminar  flame 
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speed.  While  this  expression  sometimes  gives  results  that  agree 
reasonably  well  with  experiments  for  laminar  flames,  it  does  not 
work  for  calculating  properties  of  detonations.  Some  progress,  how¬ 
ever,  has  been  made  toward  developing  a  four-step  mechanism  for 
hydrogen-oxygen  mixtures  that  is  valid  for  both  flames  and  deto¬ 
nations  [10]. 

Here  we  address  the  problem  of  developing  a  minimal  model 
that  captures  the  essential  features  of  a  flame,  a  detonation,  and 
the  deflagration-to-detonation  transition  (DDT)  in  methane-air 
mixtures.  In  previous  work,  we  described  similar  single-step  mod¬ 
els  for  low-pressure  acetylene,  low-pressure  ethylene,  and  atmo¬ 
spheric-pressure  hydrogen-air  mixtures  [11-17].  Using  these 
models  in  calculations  taught  us  a  great  deal  about  basic  physical 
processes  governing  turbulent  flame  acceleration  and  DDT,  as  sum¬ 
marized  in  [15],  and  in  high-speed  turbulence  [18].  Now  we  are 
interested  in  extending  this  approach  to  methane  and  natural  gas. 

The  work  in  this  paper  is  focused  on  developing  and  testing  a 
chemical-diffusion  model,  which,  when  coupled  to  an  appropriate 
model  for  the  fluid  dynamics,  will  be  accurate  enough  to  simulate 
the  transition  of  a  low-speed  flame  to  a  detonation  wave  in  a  large, 
confined  area  containing  a  methane-air  mixture.  To  do  this,  we 
have  taken  the  same  approach  that  was  used  previously  for  acety¬ 
lene,  ethylene,  and  hydrogen:  we  fit  parameters  so  that  they  repro¬ 
duce  experimental  and  theoretical  length  and  time  scales  of 
laminar  flames  and  detonations.  In  the  remainder  of  this  paper, 
we  describe  and  analyze  the  calibration  process  and  then  compare 
calculations  of  DDT  in  an  obstructed  channel  to  data  from  previ¬ 
ously  reported  experiments  [19,20]. 

2.  Model 


where  p,  T,  u,  v ,  P,  e,  q ,  and  Y  are  the  density,  temperature,  stream- 
wise  velocity,  transverse  velocity,  pressure,  specific  energy,  heat  re¬ 
lease,  and  fuel  mass  fraction  of  the  gas  mixture,  respectively.  The 
transport  coefficients,  viscosity  v,  mass  diffusivity  D,  and  thermal 
diffusivity  K  =  I<lp  cp,  where  I<  is  the  thermal  conductivity,  vary  with 
temperature  according  to 


t-h 

D  =  D0j,  (9) 

'j-n 

K  =  K0  — . 

P 

The  parameters  v0,  D0,  and  k0  are  assumed  to  be  constant,  and  n  is  a 
constant  chosen  to  be  0.7.  In  this  model,  the  ratio  of  specific  heats  y 
does  not  vary  with  temperature. 

We  use  a  one-step  Arrhenius  chemistry  model  such  that  the 
reaction  rate  Q  is  given  by 

Q  =  ApY  exp(-£a/PT),  (10) 

where  A  is  the  pre-exponential  factor,  R  is  the  gas  constant,  and  Ea  is 
the  activation  energy  for  the  reaction.  This  equation  is  similar  to  Eq. 
(1 )  with  a  =  1  and  b  =  0.  This  type  of  reaction  model  has  been  used  in 
past  work  for  acetylene,  ethylene,  and  hydrogen  to  solve  a  variety  of 
combustion  and  detonation  problems  involving  shock-flame  inter¬ 
actions  and  to  compute  the  properties  of  the  cellular  structure  of 
detonations  [11-14,17,16,21,15,22-24].  In  this  model,  there  is  no 
explicit  description  of  radiative  diffusion;  however,  some  radiation 
effects  could  be  implicitly  contained  in  the  adjustable  parameters. 
Energy  losses  through  boundaries  are  neglected. 


The  reactants  are  assumed  to  be  fully  premixed  and  behave  as 
an  ideal  gas,  so  that  the  flow  is  governed  by  the  compressible  reac¬ 
tive  Navier-Stokes  equations, 
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P  =  pRT, 

P  1,2  2  \ 
Pe=y-j-h^p(U2  +  V2), 


(7) 

(8) 


3.  Model  parameter  calibration 

The  one-step  Arrhenius  kinetics  used  in  this  model  cannot  ex¬ 
actly  reproduce  all  properties  of  laminar  flames  and  detonations 
in  methane-air  mixtures.  The  model  can,  however,  give  a  reason¬ 
able  approximation  of  the  key  length  and  time  scales  involved  at 
different  stages  of  DDT  [15].  During  the  first  stage,  an  initially  lam¬ 
inar  flame  is  accelerated  to  a  high-speed,  turbulent  deflagration 
wave.  The  second  stage  involves  the  formation  of  localized  regions 
of  elevated  temperatures  (hot  spots)  and  subsequent  detonation 
initiation  by  the  Zeldovich  gradient  mechanism  [25].  If  the  newly 
formed  detonation  survives,  it  spreads  to  the  rest  of  the  unburned 
material  during  the  third  stage. 

Flame  acceleration  (in  the  laboratory  frame)  occurs  primarily 
due  to  advection  by  the  induced  gas  flow,  which  can  be  orders  of 
magnitude  larger  than  the  laminar  and  turbulent  flame  speeds 
[17,16].  The  gas  flow  is  driven  by  thermal  expansion  of  the  com¬ 
bustion  products  and  increases  with  the  amount  of  heat  released 
by  the  flame  front,  which  is  a  function  of  the  surface  area  of  the 
flame.  The  flame-surface  area  increases  due  to  stretching  by  the 
flow  and  wrinkling  caused  by  turbulent  motions  and  fluid-dynamic 
instabilities.  When  the  flow  speed  approaches  the  speed  of  sound, 
shocks  form,  and  shock-flame  interactions  become  an  important 
mechanism  for  the  flame  wrinkling  and  turbulence  generation. 

Based  on  this  progression  of  physical  processes,  the  properties 
of  the  fuel-air  mixture  that  are  important  during  this  period  of 
flame  acceleration  are  the  laminar  flame  speed,  the  adiabatic  flame 
temperature,  the  viscosity,  and  the  speed  of  sound.  We  restrict  the 
discussion  to  mixtures  of  ideal  gases,  so  that  the  sound  speed  is  gi¬ 
ven  by  the  expression  c  =  ^ yRT ,  and  the  specific  heat  capacity  at 
constant  pressure  is  related  to  the  gas  constant  by  cp  =  yR/{y  -  1). 
The  temperature  rise  in  an  adiabatic  system  due  to  chemical-en¬ 
ergy  release  is  Tb  -  T0  =  qlcp  =  q(y  -  1  )/y  R.  Thus,  the  adiabatic 
flame  temperature  Th  depends  on  the  initial  system  temperature 
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T0,  the  heat  release  q ,  and  the  ratio  of  specific  heats  y.  The  laminar 
flame  speed  Si  depends  on  A,  Ea,  q,  and  k. 

Hot  spots,  which  can  evolve  to  generate  flames,  shocks,  and  det¬ 
onations,  may  arise  in  unreacted  material  from  many  types  of 
interactions  in  a  turbulent,  shock-laden  reactive  flow.  In  previous 
simulations,  we  observed  them,  for  example,  in  shock  reflections, 
vortices  behind  Mach  stems,  and  multiply-shocked  regions  of 
unreacted  material.  In  these  hot  regions,  a  spatial  temperature  gra¬ 
dient  exists,  and  the  temperature  can  be  high  enough  to  ignite  the 
reactants.  In  the  presence  of  the  temperature  gradient,  ignition  oc¬ 
curs  consecutively  in  multiple  layers  of  material  heated  to  different 
temperatures,  thus  forming  a  reaction  wave  [25].  This  wave  can 
generate  a  strong  shock  and  eventually  a  detonation.  The  survival 
of  the  newly  formed  detonation  wave  then  depends  on  local  ther¬ 
mal  and  chemical  conditions  and  geometrical  constraints  [15]. 

These  stages  of  DDT  are  controlled  by  induction  delays  behind 
strong  shocks  and  properties  of  detonation  waves.  Since  induction 
delays  correlate  with  detonation  cell  sizes,  a  model  calibrated  on 
detonation  properties  should  approximately  describe  key  phenom¬ 
ena  responsible  for  these  final  stages  of  DDT.  The  detonation  prop¬ 
erties  used  for  the  calibration  are  the  theoretical  Chapman-Jouget 
(CJ)  detonation  velocity  DCj  and  the  detonation  cell  size  2.  The 
velocity  DCJ  depends  on  the  heat  release  q  and  y.  If  the  sound  speed 
is  correct,  DQ  gives  an  indication  of  the  correct  value  of  q,  as  Tb  does 
for  a  deflagration.  Then  2  is  controlled  by  A  and  Ea. 

The  model  calibration  is  started  by  performing  a  series  of  one¬ 
dimensional  calculations  of  flame  and  detonation  structures  for  a 
range  of  input  parameters.  Then  we  choose  a  set  of  parameters  that 
most  closely  reproduces  both  the  laminar  flame  properties  and  the 
detonation  properties.  After  this,  two-dimensional  simulations  are 
used  to  compute  detonation  cell  sizes  using  the  model  parameters 
arrived  upon  in  the  previous  step.  Finally,  full  simulations  using 
the  model  parameters  are  used  to  compute  flame  acceleration 
and  DDT,  and  these  results  are  compared  to  experiments. 


3.1.  Laminar  flames 


First,  the  values  of  v0  and  k0  are  calculated  directly  from  known 
viscosity  and  thermal  conductivity  of  air  at  T=  298  K  and  P  =  1  atm. 
For  simplicity,  we  assume  the  Lewis  number  of  the  mixture  is 
equal  to  unity,  implying  D0  =  k0,  which  is  an  acceptable  approxi¬ 
mation  for  premixed  methane  and  air.  We  then  compute  proper¬ 
ties  of  a  one-dimensional  laminar  flame  by  solving  an  ordinary 
differential  equation  describing  thermal  conduction  and  energy  re¬ 
lease  inside  a  steady-state  reaction  wave, 


=  p  uc, 


dFt 

dx 


Ft  =  K7T’ 
dx 


dr 

^dx 


(11) 

(12) 


where  U  =  Sipolp  is  the  flow  velocity  in  the  frame  moving  with  the 
flame.  According  to  Eq.  (9),  the  thermal  conductivity  K  is  a  function 
of  temperature,  K=k0T°JCp,  where  Cp  =  Ry/(y-  1).  Details  of  the 
solution  method  are  given  in  Appendix  A  and  [21].  We  search  for 
a  set  of  parameters,  A,  q,  Ea,  and  y,  for  which  the  computed  value 
of  Si  matches  experimental  data  [26,27],  and  that  of  Tb  matches  val¬ 
ues  calculated  using  a  complex  reaction  mechanism  for  methane- 
air  combustion  [28].  The  laminar  flame  thickness,  X/,  is  based  on 
the  temperature  gradient,  i.e.,  X/=  (Tb  -  T0)/max|<9T/<9x|. 


3.2.  One-dimensional  detonations 

We  use  a  Zeldovich-von  Neumann-Doering  (ZND)  model  to 
compute  the  half-reaction  thickness  of  one-dimensional  detona¬ 
tions  using  the  same  reaction  model  described  in  §2  for  a  set  of 


parameters  A,  Ea,  y,  and  q.  The  reaction  zone  of  a  one-dimensional 
detonation  is  described  by 
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+  qQ , 
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where  t  is  time,  x  the  distance  from  the  shock,  U  =  Dqpolp  the  flow 
velocity  in  the  shock  frame,  c  =  y/yP/p  the  sound  speed,  and 


Dcj  =  Co 


q  p0(y2-i) 

Po  2y 


Cq  Po(y2-l A 

Po  2y  J 


(16) 


is  the  Chapman-Jouget  detonation  velocity.  The  solution  procedure 
for  this  set  of  equations  is  described  in  Appendix  B  and  [21]. 

Solutions  of  Eqs.  (13)-(16)  give  the  reaction-zone  profile,  from 
which  we  find  the  half-reaction  thickness  of  the  one-dimensional 
detonation  wave,  xd.  This  quantity  correlates  with  the  detonation 
cell  size  2  [29],  which  is  often  measured  in  detonation  experi¬ 
ments.  Even  though  detonation  cells  do  not  develop  until  after 
DDT  occurs,  k  is  generally  used  instead  of  xd  in  empirical  correla¬ 
tions  related  to  the  detonation  initiation  and  DDT  [30-41].  For 
example,  experimental  evidence  suggests  that  the  ratio  of  the  sys¬ 
tem  size  to  2  must  be  greater  than  one  for  a  sustained  detonation  to 
occur  [30].  As  we  did  for  the  laminar  flames,  we  use  measured  val¬ 
ues  of  2  [19]  (and,  equivalently,  xd)  and  values  of  DQ  calculated 
using  a  thermodynamic  equilibrium  code  [42]  to  find  a  set  of  con¬ 
sistent  model  parameters. 


3.3.  Composite  model 


The  choice  of  parameters  y  and  q  determines  Tb  for  laminar 
flames  and  DCj  for  one-dimensional  detonations.  Curves  represent¬ 
ing  values  of  q  that  give  Tb  =  2210  K  and  DCj  =  1820  m/s  as  a  func¬ 
tion  of  y  are  shown  in  Fig.  la.  We  choose  the  values  of  q  and  y  at 
the  intersection  of  these  two  curves  to  use  in  our  model.  The 
remaining  parameters,  Ea  and  A,  determine  Si  and  xd  for  fixed  k0 , 
y  and  q.  Fig.  lb  shows  the  values  of  A  for  which  5/ =  38.02  cm/s 
and  xd  =  0.229  cm  as  a  function  of  Ea.  Again,  we  choose  A  and  Ea 
at  the  point  of  intersection  of  these  curves.  The  results  of  this 
one-dimensional  calibration  process  are  given  in  Table  1  for  a  stoi¬ 
chiometric  (9.5%)  methane-air  mixture.  The  computed  input 
parameters  and  the  resulting  output  for  Sb  Tb,  DQ,  and  2  are  shown 
along  with  values  of  these  parameters  from  the  literature.  We  note 
that  the  determination  of  2  for  high-activation-energy  mixtures  is 
not  precise  because  the  detonation  cells  are  highly  irregular. 
Hence,  ranges  for  2  and  xd  are  given.  The  value  of  xd  =  0.229  cm 
used  in  the  calibration  process  was  chosen  as  a  representative  va¬ 
lue  within  the  range  shown  in  Table  1.  In  practice,  a  range  of  Ea  and 
A  would  give  satisfactory  values  of  2. 


3.4.  Two-dimensional  detonation  cells 


High-resolution,  two-dimensional  simulations  of  a  detonation 
propagating  in  a  planar  channel  [43]  were  performed  using  the 
model  parameters  chosen  in  the  previous  section.  The  computa¬ 
tional  cell  size  at  the  detonation  front  was  chosen  so  that  234  cells 
spanned  xd.  This  level  of  grid  refinement  was  necessary  to  properly 
resolve  the  transverse  wave  structures  that  generate  the  character¬ 
istic  cellular  pattern  as  an  unstable  detonation  propagates.  The 
detonation  cells  that  formed  were  highly  irregular,  ranging  in  size 
from  10  to  20  cm.  Fig.  2  shows  a  sequence  of  density  contours  over 
a  distance  equivalent  to  several  detonation  cell  widths.  In  this 
example,  the  maximum  distance  between  two  adjacent  triple 
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Fig.  1.  Parametric  curves  for  which  (a)  Tb  =  2210  K,  DCj  =  1820  m/s  and  (b)  S/  =  38.02 
RT0,  y,  EJRTq,  and  A  used  in  conjunction  with  the  reaction  model  (Eq.  (10)). 


(b) 

Xd  =  0.229  cm.  The  points  of  intersection  in  the  two  figures  give  the  values  of  qM/ 


Table  1 

Input  model  parameters  and  computed  properties  of  reaction  waves  for  stoichiom¬ 
etric  methane-air  mixture. 


Input 

Po 

1  atm 

To 

298  K 

M 

27  g/mol 

y 

1.197 

A 

1.64  x  1013  cm3/g  s 

Ea 

67.55 RT0 

Q 

39.0  RTo/M 

Vo 

3.6  x  10-6  g/s  cm  K0-7 

Ko  =  Do 

6.25  x  10-6  g/s  cm  K0  7 

Output 

Calculated  values 

Target  values 

s, 

38.02  cm/s 

34-45  cm/s  [26,27] 

Tb 

2210  K 

2200-2230 K [28] 

xf 

0.0439  cm 

Dcj 

1820  m/s 

-1815  m/s  [42] 

xd  (2) 

0.229  cm  (16-23  cm)* 

0.13-0.62  cm  *  (13-31  cm)  [19] 

Based  on  estimated  ratio  of  50  <  2/xd  <100  [29]. 


points  is  approximately  16  cm,  but  this  distance  varies  as  the  det¬ 
onation  propagates  and  new  triple  points  form,  giving  rise  to  the 


range  of  cell  sizes  discussed  above.  The  sizes  of  the  computed  cells 
are  consistent  with  the  range  13-31  cm  observed  in  experiments 
[19].  The  computed  ratio  of  2/xd  for  this  model  is  in  the  range  of 
43-87,  which  is  consistent  with  previous  calculations  [29].  Thus, 
the  chemistry  model  with  coefficients  calibrated  using  one-dimen¬ 
sional  models  produces  two-dimensional  detonations  with  length 
and  time  scales  consistent  with  experimental  observations. 

4.  Two-dimensional  channels  with  obstacles 

We  next  use  this  calibrated  model  to  calculate  multidimen¬ 
sional,  turbulent,  accelerating  flames  and  subsequent  DDT.  Ob¬ 
structed  channels  promote  faster  flame  acceleration  than  smooth 
channels  and,  hence,  provide  a  convenient  testbed  for  studying 
DDT.  In  addition,  this  configuration  has  been  used  for  several 
experimental  studies,  and  so  there  is  some  data  that  can  be  used 
for  comparisons. 

The  model  planar  2D  channel  used  in  the  calculations  is  shown 
in  Fig.  3.  The  channel  is  closed  at  the  left  end  (x  =  0)  and  can  either 
be  open  to  the  atmosphere  or  closed  at  the  opposite  end  (x  =  L).  At 
the  bottom  plane  (y  =  0)  and  the  x  =  0  boundary,  we  assume  a  non¬ 
slip,  adiabatic  wall,  i.e.,  u  =  v  =  0  and  dT/dn  =  dY/dn  =  0,  where  n  de¬ 
notes  the  direction  normal  to  the  surface.  The  same  nonslip 


■  m  i  _ Li. _ mm ,  „i  _ umm tu _ _ ^mm _ :jm  j_. ^ 

450  460  460  470  470  480  480  480  490  490 

x  (cm)  x  (cm)  x  (cm)  x  (cm)  x  (cm)  x  (cm)  x  (cm) 


Density  (g/cm3) 

0.0005  0.0015  0.0025  0.0035  0.0045  0.0055  0.0065  0.007 

Fig.  2.  Selected  density  maps  near  the  reaction  front  of  a  detonation  propagating  through  a  32  cm  wide  by  1024  cm  long  channel  filled  with  a  stoichiometric  mixture  of 
methane  and  air  calculated  using  the  reaction  model  (Eq.  (10))  and  the  calibrated  parameters  listed  in  Table  1  at  times  t  =  (a)  2.480  ms,  (b)  2.505  ms,  (c)  2.542  ms,  (d) 
2.586  ms,  (e)  2.620  ms,  (f)  2.648  ms,  and  (g)  2.677  ms. 
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Fig.  3.  Computational  setup.  Obstacles  are  evenly  spaced  along  the  entire  length  of 
the  channel.  Walls  and  obstacle  surfaces  are  adiabatic  no-slip  reflecting  boundaries. 
Initial  flame  radius  is  0.25  cm. 


boundary  conditions  are  imposed  along  the  face  of  every  obstacle 
and  at  the  right  boundary  x  =  L  for  a  closed  channel.  For  an  open 
channel,  a  zero-gradient  outflow  boundary  condition  was  imposed 
at  x  =  L.  We  assume  the  channel  is  symmetric  and  simulate  only 
the  lower  half,  so  symmetry  conditions  du/dy  =  dY/d  y  =  d T/ 
dy  =  dP/dy  =  v  =  0  are  applied  at  the  channel  center  line,  y  =  d/2. 
The  obstacles  are  taken  to  be  2  cm  thick,  and  their  heights  are 
set  based  on  the  desired  blockage  ratio,  £  =  2 fr/d.  Obstacle  spacings 
S  are  set  equal  to  d. 

We  consider  three  different  configurations  chosen  to  be  similar 
to  experimental  systems  [19,20]  and  summarized  in  Table  2.  In  the 
first  configuration,  L  =  216.2  cm  and  d  =  7.6  cm,  which  models  the 
7.6  x  7.6  cm  square  channel  used  in  [20].  In  those  experiments, 
both  ends  of  the  channel  were  closed,  and  the  initial  pressure  in 
the  unburned  gas  mixture  was  47  kPa.  The  second  is  a  slightly  lar¬ 
ger  channel  with  d  =  1 7.4  cm  and  L  =  1187.8  cm,  which  is  similar  to 
the  circular  cross-section  tube  (diameter  17.4  cm)  used  in  [19]. 
Here,  the  right  end  of  the  channel  is  open  to  the  atmosphere,  and 
the  initial  gas  pressure  is  atmospheric,  as  in  the  experiments.  The 
third  configuration  (d  =  52cm,  L  =  2130  cm)  is  similar  to  the 
52  cm  diameter  tube  used  in  [19],  where  the  channel  is  also  open 
to  the  atmosphere  at  x  =  L.  For  each  test  case,  the  channel  is  uni¬ 
formly  filled  with  a  stoichiometric  methane-air  mixture. 

To  ignite  the  mixture,  we  place  a  quarter-circular  region  of  hot, 
burned  material  at  the  left  wall  on  the  centerline  and  add  a  small 
amount  of  extra  energy  to  the  burned  region.  The  additional  energy 
per  unit  mass  is  on  the  order  of  the  chemical-energy  release  q , 
which  could  model  ignition  by  a  low-energy  (~100mJ)  spark. 
The  resulting  weak  shock  wave  is  not  nearly  strong  enough  to 
ignite  a  detonation  directly.  It  only  causes  multiple  shock  reflec¬ 
tions  and  shock-flame  interactions  that  distort  and  wrinkle  the 
flame  front. 

Eqs.  (2)-(8)  are  solved  using  an  explicit,  second-order,  Godu¬ 
nov-type  numerical  scheme  incorporating  a  Riemann  solver.  The 
integration  is  performed  on  a  structured  adaptive  mesh  based  on 
the  fully  threaded  tree  data  structure  [44].  The  mesh  refinement 
is  dynamically  controlled  by  gradients  of  density,  temperature, 
and  composition.  Typically,  the  maximum  computational  cell  size 
(away  from  shocks  and  flame  fronts)  is  dxmin  =  0.29  cm,  and  the 
minimum  refined  cell  size,  dxmin  =  0.0163  cm,  which  corresponds 
to  3-4  computational  cells  per  laminar  flame  thickness. 

Fig.  4  shows  the  one-dimensional  flame  structure  calculated 
from  a  two-dimensional  simulation  of  a  planar  flame  propagating 


Table  2 

Model  configurations. 


Configuration* 

7.6 

17.4 

52 

L  (cm) 

216.2 

1187.8 

2130 

d  (cm) 

7.6 

17.4 

52 

i 

1/3,  2/3 

0.3,  0.6 

0.3,  0.6 

x  =  l  boundary 

Closed 

Open 

Open 

Po(atm) 

0.464 

1 

1 

Experimental  channel  cross-section 

Square 

Circular 

Circular 

Experiments 

[20] 

[19] 

[19] 

Note  that  tube  diameter  (for  circular  cross-sections)  or  side  length  (for  square 
cross-sections)  is  used  as  a  naming  convention  for  each  experiment. 


Fig.  4.  Temperature  (solid)  and  reaction-rate  (dashed)  profiles  calculated  using  Eqs. 
(11),  (12)  (lines),  two-dimensional  Navier-Stokes  equations  with  dxmin  = 
0.018125  cm  (triangles),  and  the  two-dimensional  Navier-Stokes  equations  with 
dxmin  =  0.001 13  cm  (squares). 


in  a  smooth  channel  for  this  resolution  (large  symbols)  and  a  finer 
resolution,  dxmin  =  0.00113  cm  (small  symbols).  The  flame  struc¬ 
ture  calculated  using  the  steady-state  one-dimensional  laminar 
flame  model  (Eqs.  (11),  (12))  is  also  shown  in  Fig.  4  (lines).  The 
high-resolution  two-dimensional  calculation  reproduces  the  theo¬ 
retical  flame  structure  and  gives  nearly  the  same  flame  speed. 
There  are  some  differences  in  the  reaction-rate  profiles  calculated 
at  the  lower  resolution  ( dxmin  =  0.0163  cm).  The  computed  laminar 
flame  speed  is  approximately  12%  smaller  than  the  theoretical  lam¬ 
inar  flame  speed.  More  resolution  tests  are  discussed  in  detail  in 
Section  4.4.  Here  we  only  note  that  in  larger  channels,  dxmin  is  lim¬ 
ited  by  the  available  computational  resources,  and  in  most  cases, 
the  simulations  must  be  somewhat  under-resolved. 


4A.  Configuration  7.6 

Johansen  and  Ciccarelli  [20]  examined  the  development  and 
acceleration  of  a  turbulent  flame  in  a  7.6  x  7.6  cm  square  cross- 
section  channel.  Obstacles,  of  heights  1.27, 1.9,  and  2.53  cm  corre¬ 
sponding  to  f  =  1/3, 1/2,  and  2/3,  respectively,  were  spaced  7.6  cm 
apart  over  the  entire  244  cm  channel  on  both  the  top  and  bottom 
walls.  Both  ends  of  the  chamber  were  closed,  and  the  pressure  of 
the  stoichiometric  methane-air  mixture  inside  the  chamber  was 
initially  47  kPa.  The  mixture  was  ignited  by  an  electric  spark  at 
the  centerline  of  the  channel. 

We  have  simulated  a  two-dimensional  rectangular  channel 
(Fig.  3)  with  obstacle  spacings  and  heights  identical  to  those  in 
the  experiments.  After  ignition,  we  track  the  position  and  velocity 
of  the  leading  edge  of  the  flame  front  as  well  as  the  total  length  of 
flame  surface  created  as  the  flame  evolves.  Flame  velocities  are 
computed  at  discrete  locations  along  the  length  of  the  channel 
and  represent  an  average  velocity  over  the  interval  between  two 
successive  measurement  locations.  The  flame  surface  is  calculated 
by  summing  the  total  length  of  the  isosurface  on  which  Y  =  0.5  at  a 
particular  instance  in  time.  Fig.  5  compares  measured  [20]  and 
computed  flame  velocities  and  flame  surface  areas  for  £  =  1  /3  and 
£  =  2/3.  At  early  times  in  the  flame  development,  the  simulations 
and  experiments  show  similar  flame  velocities  for  both  blockage 
ratios.  Differences  arise  further  downstream  (x~  100-150  cm)  as 
the  flame  evolves. 

The  flame  acceleration  process  occurs  in  three  phases,  each 
of  which  can  be  characterized  by  the  dominant  mechanism 
driving  the  growth  in  flame  surface  area.  In  the  first  phase,  the 
flame  is  folded  and  stretched  by  a  laminar  flow  field  that  is 
induced  by  the  thermal  expansion  of  the  combustion  products. 
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(a)  (b) 


Fig.  5.  Configuration  7.6:  (a  and  b)  Computed  and  measured  [20]  flame-propagation  velocities  and  (c)  computed  flame-surface  as  a  function  of  the  position  of  the  leading 
edge  of  the  reaction  front. 


The  development  of  flame  surface  is  similar  for  £  =  1/3  and  2/3  for 
this  phase,  x  <  70  cm  (Fig.  5c).  Temperature  maps  of  the  leading 
edge  of  the  flame  as  it  passes  over  the  first  obstacle  are  shown  in 
Fig.  6a  and  b.  There  is  relatively  little  wrinkling  of  the  flame  front 
during  this  phase. 

In  the  second  phase,  the  predominant  mechanism  for  increasing 
the  total  length  of  flame  surface  is  wrinkling  by  fluid  dynamic 
instabilities  (e.g.,  Kelvin-Helmholtz  and  Rayleigh-Taylor)  and  tur¬ 
bulent  fluctuations.  Localized  regions  of  vorticity  stretch  and  frag¬ 
ment  a  continuous  flame  front,  thereby  increasing  the  total 
amount  of  flame-surface  area.  The  energy  released  at  the  flame 
surface  causes  the  thermal  expansion  of  the  product  gases,  which 
causes  a  net  flow  through  the  channel.  Shear  layers  develop  down¬ 
stream  of  obstacles  as  the  fluid  is  accelerated  through  the  re¬ 
stricted  cross-sectional  area  above  them.  Fluid  dynamic 
instabilities  and  turbulence  in  the  shear  layers  contribute  to  gener¬ 
ation  of  more  flame-surface  area.  The  flame  surface  and  velocity 
become  substantially  larger  for  the  £  =  2/3  case  than  for  the  {  =  1/ 
3  case  for  x  >  70  cm.  In  the  £  =  2/3  case,  the  shear  layers  develop 
more  quickly,  since  the  flow  is  accelerated  to  a  higher  velocity  in 
the  smaller  gap  between  obstacles. 

It  may  be  possible  to  describe  the  maximum  turbulent  flame 
speed  attained  during  this  stage  of  flame  acceleration  using  a 
one-dimensional  model,  such  as  that  described  by  Bradley  et  al. 
[45,46];  however,  the  turbulent  flame  brushes  developed  in  our 
calculations  are  qualitatively  different  from  the  idealized  one¬ 
dimensional  turbulent  flame.  The  extended  burning  regions  be¬ 


tween  obstacles  far  behind  the  leading  edge  of  the  flame  release 
significant  amounts  of  energy,  and  in  some  regions  the  flame  prop¬ 
agates  normal  to  the  direction  of  the  induced  gas  flow. 

Temperature  maps  of  the  flame  fronts  near  x  =  100  cm  (Fig.  6c 
and  d)  show  that  there  is  much  more  small-scale  flame  structure 
and  therefore  much  more  flame  surface  in  the  £  =  2/3  case.  The  in¬ 
crease  in  total  flame-surface  area  (see  Fig.  5c)  continues  as  long  as 
substantial  amounts  of  fuel  remain  in  between  obstacles  behind 
the  foremost  part  of  the  flame  front.  Then,  when  most  of  this  fuel 
is  depleted,  the  amount  of  flame-surface  decreases  rapidly.  The 
maximum  flame-surface  area  developed  for  £  =  2/3  is  larger  than 
that  for  {  =  1  / 3  because  the  increased  velocity  of  the  flame  front  al¬ 
lows  the  flame  to  propagate  farther  in  the  channel  before  these  ex¬ 
tended  reaction  zones  can  burn  out.  The  extra  amount  of  flame 
surface  present  in  these  regions  results  in  faster  depletion  of  the 
fuel,  and  hence  the  steep  drop  in  flame-surface  area  for 
x  >  130  cm. 

This  rapid  decline  in  flame-surface  area  is  slowed  as  the  flame 
enters  the  third  phase.  When  the  speed  of  the  induced  flow  ap¬ 
proaches  the  speed  of  sound  in  the  unburned  mixture  ahead  of 
the  flame,  energy  released  at  the  leading  edge  of  the  flame  front 
generates  weak  pressure  waves  that  propagate  ahead  of  the  flame. 
These  pressure  waves  later  become  shocks,  which  can  reflect  from 
obstacles  and  walls.  The  reflected  shocks  also  collide  and  interact 
with  portions  of  the  reaction  front.  The  flame  surface  is  wrinkled 
by  these  shock-flame  interactions  that  promote  Richtmyer-Mesh- 
kov  instabilities.  The  turbulence  generated  by  these  instabilities  is 
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Fig.  6.  Temperature  maps  near  the  leading  edge  of  the  flame  for  configuration  7.6  with  £  =  1/3  (left)  and  £  =  2/3  (right)  at  several  locations  throughout  the  channel:  (a  and  b) 
near  first  obstacle,  (c  and  d)  near  x  =  100  cm,  and  (e  and  f)  near  x  =  200  cm.  The  top  temperature  scale  is  for  burned  material,  and  the  bottom  scale  is  for  unburned  material. 
Time  increases  from  top  to  bottom  in  each  column. 


not,  necessarily,  homogeneous,  isotropic,  Kolmogorov  turbulence. 
The  nature  of  this  turbulence  and  its  interaction  with  the  flame 
are  interesting  areas  of  future  investigation.  Here,  the  additional 
energy  release  caused  by  flame-surface  wrinkling  helps  to  sustain 
the  flame  speed  and  slow  the  decline  in  the  net  flame  surface.  This 
process  occurs  for  x  >  1 60  cm  for  £  =  2/3,  but  does  not  begin  to  oc¬ 
cur  in  the  £  =  1  /3  case  until  the  flame  reaches  the  end  of  the  do¬ 
main.  A  longer  channel  would  be  necessary  to  observe  significant 
shock-flame  interactions  for  the  £  =  1  /3  case.  Fig.  6e  and  f  compare 
the  temperature  maps  for  £  =  1/3  and  2/3  when  the  leading  edge  of 
the  flame  is  near  x  =  200  cm.  A  well-defined  shock  wave  has 
formed  in  the  £  =  2/3  case,  but  the  waves  ahead  of  the  flame  front 
have  not  yet  coalesced  into  a  shock  for  £  =  1/3. 

4.2.  Configurations  1 7.4  and  52 

Kuznetsov  et  al.  [19]  performed  similar  experiments  in  circular 
cross-section  tubes  with  diameters  of  17.4  cm  and  52.0  cm.  The 
obstacles  in  the  tubes  were  annular  orifice  plates  that  were  spaced 
one  diameter  apart.  The  blockage  ratio  defined  in  these  experi¬ 
ments  is  then  =  1  —  (D*/D)2,  where  D  and  D*  are  the  tube  diame¬ 
ter  and  orifice  diameter,  respectively.  In  these  experiments,  one 
end  of  the  tube  was  left  open  to  the  atmosphere,  and  the  initial 
gas  pressure  was  atmospheric  throughout  the  tube.  They  ignited 
the  uniform  stoichiometric  methane-air  mixture  near  the  tube  axis 
at  the  closed  end.  Photodiodes  were  placed  at  various  positions 
along  the  walls  of  the  tubes,  and  reaction-front  velocities  were  cal¬ 
culated  based  on  time-of-arrival  measurements. 

After  the  initial  acceleration  period,  two  propagation  velocity 
regimes  were  found.  The  first,  commonly  referred  to  as  the  “chok¬ 
ing”  regime,  is  characterized  by  a  velocity  close  to  1/2DC;  [47].  The 
second  regime,  the  “quasi-detonation,”  is  characterized  by  a  flame- 
front  velocity  just  less  than  DCj.  For  a  blockage  ratio  of  <f  =  0.3,  the 
experimentally  measured  flame  speed  fluctuates  between  the 
speeds  typical  of  these  two  propagation  regimes.  This  indicates 
that  the  17.4  cm  diameter  tube  is  close  to  the  critical  size  for  det¬ 
onation  propagation,  as  supported  by  the  observation  that  D/2  <  1 
for  the  stoichiometric  methane-air  system.  For  the  larger  blockage 
ratio,  <jf  =  0.6,  the  experimental  flame  velocity  approaches  a  steady 
value  of  approximately  700  m/s,  a  velocity  characteristic  of  the 


choking  propagation  regime.  Similar  results  were  obtained  for 
the  D  =  52.0  cm  tube. 

We  performed  simulations  similar  to  these  experiments  using 
the  configuration  shown  in  Fig.  3  with  d  =  D  and  £  =  <f  for  both 
the  17.4  cm  (configuration  17.4)  and  52  cm  (configuration  52) 
cases.  Because  of  the  differences  in  geometry,  the  obstacle  heights 
in  the  simulations,  h ,  are  slightly  larger  than  the  heights  of  the  ori¬ 
fice  plates  in  the  experiments,  \i  =  (D  -  D*)/2,  for  the  same  block¬ 
age  ratios.  Fig.  7  shows  the  calculated  flame  velocities  and 
surface  areas  for  configuration  17.4  as  a  function  of  the  position 
of  the  leading  edge  of  the  reaction  wave  for  £  =  0.3  and  0.6.  The 
symbols  on  Fig.  7a, b  are  flame-velocity  data  from  the  experiments 
in  the  17.4  cm  tubes  [19].  In  both  the  calculations  and  experi¬ 
ments,  the  flame  accelerates  to  a  fixed  velocity  characteristic  of 
the  choking  regime,  which  then  either  undergoes  DDT  or  continues 
to  propagate  at  this  average  speed.  The  initial  flame  acceleration  is 
similar  to  that  described  in  the  previous  section:  flames  are 
stretched  by  the  thermal-expansion-induced  flow,  wrinkled  and 
torn  by  turbulence  and  fluid  dynamic  instabilities,  and  further  frag¬ 
mented  by  shock-flame  interactions.  For  example,  the  black  line  in 
each  frame  in  Fig.  8  shows  the  progress  of  a  shock  colliding  with 
and  passing  through  a  flame.  As  the  shock  passes  through,  signifi¬ 
cantly  more  flame-surface  area  is  created  behind  it.  For  the  cases 
shown  in  Fig.  7,  the  channel  is  long  enough  for  the  flames  to  pro¬ 
gress  through  all  three  stages  of  the  acceleration  process. 

The  evolution  of  flame  surfaces  for  0  <  x  <  450  cm  shown  in 
Fig.  7c  follows  the  same  trend  as  in  the  early  stages  of  configura¬ 
tion  7.6.  For  x  >  450  cm  and  £  =  0.3,  the  flame-surface  area  shar¬ 
ply  decreases.  At  this  point,  however,  the  reaction-front  velocity 
jumps  to  DCj ,  indicating  that  a  detonation  was  initiated  and  sur¬ 
vived.  The  sequence  of  events  that  lead  to  DDT  is  shown  in 
Fig.  9a-e.  Strong  shock  waves  formed  ahead  of  the  flame  front  re¬ 
flect  off  of  the  channel  wall  and  the  faces  of  obstacles,  which  re¬ 
sults  in  even  stronger  waves  and  more  shock-flame  interactions. 
Eventually,  Mach  stems  form  and  these  raise  the  local  temperature 
close  to  the  ignition  point.  These  regions  of  elevated  temperature, 
or  hot  spots,  may  or  may  not  ignite,  depending  on  the  ignition  de¬ 
lay  time  of  the  mixture  and  the  length  of  time  the  temperature  re¬ 
mains  elevated.  The  hot  spot  created  by  a  Mach  stem  just 
beginning  to  reflect  from  the  base  of  the  obstacle  in  Fig.  9b  (at 
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Fig.  7.  Configuration  17.4:  (a  and  b)  Reaction  front-propagation  velocities  and  (c)  calculated  flame-surface  lengths  as  a  function  of  the  position  of  the  leading  edge  of  the 
reaction  front. 


450  cm)  leads  to  a  successful  detonation  ignition.  The  detonation 
then  propagates  into  unburned  fuel,  catches  up  to  the  leading 
shock  wave,  and  proceeds  to  consume  nearly  all  unburnt  fuel  in 
that  region  (Fig.  9c-e).  This  is  qualitatively  similar  to  the  process 
reported  for  detonation  ignition  in  hydrogen-oxygen  mixtures  in 
channels  with  obstacles  [16,17].  The  size  of  the  system,  however, 
is  considerably  larger  for  the  methane-air  mixture. 

The  velocity  curves  shown  in  Fig.  7a  for  d  =  17.4  cm  and  {  =  0.3 
indicate  that  the  detonation  propagates  at  a  speed  much  less  than 
Dq.  This  is  a  result  of  recurring  detonation  diffractions  that  contin¬ 
ually  decouple  the  flame  from  the  leading  shock.  The  detonation  is 
then  reignited  at  a  subsequent  obstacle  in  the  manner  discussed 
above.  Although  a  detonation  propagates  at  a  speed  greater  than 
or  equal  to  DQ,  a  fast  deglagration  (decoupled  flame  and  shock) 
propagates  significantly  slower.  Thus,  the  time-averaged  velocity 
for  this  quasi-detonation,  characterized  by  intermittent  periods 
of  detonation  and  fast  deflagration  propagation,  is  below  DCj.  An 
example  of  detonation  failure  and  subsequent  reignition  is  shown 
in  Fig.  9e-j.  This  repeated  ignition  and  decoupling  process  leads  to 
the  observed  smaller  propagation  velocities  for  this  case. 

The  velocities  obtained  for  the  d  =  52  cm  cases  ({  =  0.3  and  0.6) 
are  shown  in  Fig.  10.  For  £  =  0.3,  the  computed  and  measured  [19] 
flame  velocities  are  very  similar  during  the  initial  flame  accelera¬ 
tion  period  (x  <  700  cm).  Near  700  cm,  DDT  occurs  in  the  simu¬ 
lated  system  and  the  computed  velocity  jumps  to  ~1 800  m/s. 
This  propagation  speed  is  much  closer  to  DCj  than  that  observed 
in  configuration  17.4,  £  =  0.3  since  fewer  instances  of  shock-flame 


decoupling  take  place  in  the  larger  channel.  In  the  experiments, 
DDT  first  occurs  farther  downstream,  near  x~  1000  cm,  and  the 
quasi-detonation  velocity  is  somewhat  smaller  than  the  calculated 
value.  For  £  =  0.6,  the  computed  flame  acceleration  is  close  to  the 
experimental  data.  In  the  simulations,  several  instances  of  DDT 
were  observed,  while  no  DDT  occurred  in  the  experiments.  The  to¬ 
tal  flame  surface  for  the  £  =  0.6  simulation  is  everywhere  greater 
than  that  of  the  f  =  0.3  simulation.  Larger  pockets  of  unburned  fuel 
between  obstacles  take  longer  to  burn  and  delay  the  onset  of  the 
rapid  decline  in  flame-surface  area  that  occurred  in  configurations 
7.6  and  17.4  for  the  larger  blockage  ratios.  By  the  time  the  pockets 
of  fuel  begin  to  burn  out,  the  leading  edge  of  the  propagating  flame 
has  already  accelerated  to  the  point  where  frequent  shock-flame 
interactions  significantly  increase  the  amount  of  flame  surface, 
leading  to  less  rapid  net  losses  of  flame-surface  area. 

The  detonation  is  less  likely  to  fail  when  the  orifice  diameter,  D*, 
is  large  compared  to  the  detonation  cell  size.  Peraldi  et  al.  [30]  sug¬ 
gested  that  a  suitable  criterion  for  whether  or  not  DDT  can  occur  in 
obstructed  channels  is  D*/A  >  1.  Dorofeev  et  al.  [48]  proposed  a  dif¬ 
ferent  metric  based  on  a  length  scale  that  depends  on  the  distance 
between  the  leading  edges  of  adjacent  obstacles,  5,  (cf.  Fig.  3)  and 
the  orifice  diameter.  They  showed  that  the  metric  L*/2  >  7,  where 
L*  =  (5  + D)/2/(l  -  D*/D)  was  a  reliable  indicator.  Table  3  shows 
approximate  values  for  D*// 1  and  L*/2  for  each  of  the  experimental 
systems. 

For  the  52  cm  tube  with  £  =  0.3,  D*/2  and  L*/2  are  both  much  lar¬ 
ger  than  the  critical  values.  The  experimental  data  show  that  DDT 
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Fig.  8.  Temperature  maps  near  the  leading  edge  of  flame  in  configuration  17.4  with  £  =  0.3  as  a  shock  interacts  with  the  flame  front.  The  heavy  black  line  indicates  the  location 
of  a  shock,  and  the  arrow  indicates  the  direction  of  propagation.  The  top  temperature  scale  is  for  burned  material,  and  the  bottom  scale  is  for  unburned  material.  Time 
increases  from  top  to  bottom. 
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Fig.  9.  Configuration  17.4,  £  =  0.3:  Temperature  maps  near  the  leading  edge  of  the  reaction  front  that  show  DDT  (left  column)  and  shock-flame  decoupling  followed  by 
subsequent  detonation  reignition  (right  column).  The  top  temperature  scale  is  for  burned  material,  and  the  bottom  scale  is  for  unburned  material.  Time  increases  in 
alphabetical  order  in  panels  (a)  through  (j). 
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Fig.  10.  Configuration  52:  (a  and  b)  Reaction  front-propagation  velocities  and  (c)  calculated  flame-surface  lengths  as  a  function  of  the  position  of  the  leading  edge  of  the 
reaction  front. 


Table  3 

Values  of  DDT  criterion  proposed  by  Peraldi  et  al.  [30]  and  Dorofeev  et  al.  [48]  for  the 
experimental  systems  [19]  and  the  model  configurations  used  in  the  simulations. 


£ 

Experiments 

Simulations 

D  (cm) 

DfX 

l\X 

DDT? 

d  (cm) 

D'/2 

I'M 

DDT? 

0.3 

17.4 

0.75 

5.6 

Yes 

17.4 

0.6 

3 

Yes 

52 

2.25 

17 

Yes 

52 

1.9 

9 

Yes 

0.6 

17.4 

0.6 

2.5 

No 

17.4 

0.4 

1.5 

No 

52 

1.7 

7 

No 

52 

1.1 

4.5 

Yes 

occurred  (Fig.  10a),  but  the  propagation  velocity  was  lower  than 
DCj  and  consistent  with  the  quasi-detonation  velocity  of  between 
1400  and  1500  m/s.  When  £  =  0.6  for  the  52  cm  tube,  DDT  should 
occur  based  on  the  D*/2  criterion,  but  the  L*/2  criterion  suggests 
that  the  system  is  near  the  critical  value.  The  experiments  show 
that  DDT  did  not  occur  for  this  case  (Fig.  10b). 

For  the  D  =  17.4  cm  tube,  all  of  the  criteria  are  less  than  their 
critical  values  (Table  3),  yet  there  is  some  indication  of  DDT  in 
the  experimental  data  for  £  =  0.3  (Fig.  7a).  In  this  case,  the  propaga¬ 
tion  velocity  alternates  between  quasi-detonation  and  choking, 
indicating  long  periods  of  propagation  as  a  fast  deflagration  be¬ 
tween  shorter  periods  of  detonation  propagation. 

In  the  simulations,  the  orifice  size  D =  d  -  2h  is  always  smaller 
than  D*  in  the  experiments  when  d  =  D  and  Accordingly,  D'/2 
and  L/2,  where  L  =  d/(  1  -  D  /d),  are  also  smaller  than  D*/2  and  L*/2. 
Based  on  criteria  D'/2  >  1  and  1/2  >  7  and  assuming,  for  the  mo¬ 
ment,  that  the  characteristic  2  for  the  simulations  is  the  same  as 


that  for  the  experiments  (2=19  cm),  the  only  case  that  should  be 
expected  to  undergo  DDT  is  configuration  52  with  { =  0.3. 
Fig.  10a  does  show  DDT  for  this  case,  and  the  resulting  propagation 
velocity  is  larger  than  the  experimental  data.  The  velocity  data 
from  configuration  52  with  £  =  0.6  also  show  occasional  transitions 
to  detonation,  but  the  average  propagation  speed  is  much  smaller 
than  that  obtained  in  the  case  of  smaller  blockage  ratio.  For  config¬ 
uration  17.4,  DDT  occurs  in  the  ^  =  0.3  case  in  spite  of  the  low  val¬ 
ues  of  D/2  and  L'/2. 

The  simulation  data  suggest  that  the  transition  criteria  D/2  >  1 
and  L'/2  >  7  may  not  be  sufficient  to  predict  whether  or  not  DDT 
will  occur  in  the  systems  under  consideration.  Transitions  to  deto¬ 
nations  were  found  for  D'/2  and  L'/2  as  low  as  0.6  and  3,  respec¬ 
tively;  however,  some  care  must  be  taken  when  evaluating  these 
transition  criteria.  First,  detonation  cells  in  methane-air  mixtures 
are  highly  irregular  structures,  so  the  specification  of  an  average 
2  has  large  uncertainties,  thus  making  it  difficult  to  compare  the 
sizes  of  measured  and  calculated  detonation  cells.  It  is  possible  that 
the  characteristic  size  of  the  simulated  2  is  smaller  than  that  used 
in  computing  the  ratios,  which  would  have  the  effect  of  increasing 
D  /2  and  L  /2  for  each  simulation.  Second,  the  differences  in  geom¬ 
etry  between  the  simulations  and  experiments,  particularly  the 
differences  in  obstacle  heights,  could  produce  different  shock 
reflections,  and  this  can  affect  the  occurrence,  time,  and  location 
of  DDT.  Given  these  differences,  there  is  no  reason  to  expect  that 
the  critical  values  of  the  empirical  transition  criteria  for  tubes  will 
be  the  same  as  those  for  a  two-dimensional  channel.  Any  further 
consideration  of  the  matter  would  require  a  detailed  study  of  the 
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cellular  structure  of  stoichiometric  methane-air  detonations  in  or¬ 
der  to  obtain  a  truly  reliable  estimate  of  the  average  cell  size. 

4.3.  Sensitivity  analysis 

It  is  important  to  evaluate  the  sensitivity  of  the  computed  flame 
acceleration  and  occurrences  of  DDT  to  variations  in  the  parame¬ 
ters  of  the  chemistry  model.  We  first  consider  the  sensitivity  of 
the  flame  acceleration  by  evaluating  how  much  velocity  profiles 
change  when  there  are  systematic  variations  in  the  length  and  time 
scales  of  the  laminar  flame.  This  is  done  by  changing  model  param¬ 
eters  A,  q ,  Ea,  and  y  to  create  moderate  (10-15%)  variations  in  the 
laminar  flame  velocity,  adiabatic  flame  temperature,  and  specific- 
heat  ratio.  We  next  consider  the  impact  on  DDT  of  changing  these 
parameters.  As  discussed  earlier,  model  parameters  also  affect  var¬ 
ious  length  and  time  scales  of  detonation  waves,  in  particular,  DCj, 
A,  and  ignition  delay  times. 


The  first  question  to  address  is  how  changes  in  model  parame¬ 
ters  affect  the  propagation  velocity  of  the  flame  in  an  obstructed 
channel.  We  performed  a  series  of  simulations  using  configuration 
17.4  with  {  =  0.3  and  0.6.  Fig.  11  shows  the  velocity  of  the  leading 
edge  of  the  flame  front  as  a  function  of  position.  In  Fig.  1  la,  we  con¬ 
sider  several  different  sets  of  parameters  that  yield  the  same 
Tb  =  2210 1<  but  different  values  of  the  one-dimensional  laminar 
flame  speed.  A  description  of  each  parameter  set,  is  given  in  Ta¬ 
ble  4  along  with  the  corresponding  laminar  flame  and  ID  detona¬ 
tion  properties  calculated  from  Eqs.  (11 )— (16).  The  solid  lines 
represent  calculations  for  a  parameter  set  0>\  that  gives 
5/  =  38.02  cm/s,  while  the  laminar  burning  velocity  of  the  flames 
represented  by  the  dashed  (^2)  and  dotted  (^3)  lines  is  32.5  cm/ 
s.  In  ^2,  the  decrease  in  S/  was  caused  by  reducing  A  to 
1.2  x  1013cm3/gs.  In  ^3,  the  same  reduction  in  Sj  was  brought 
about  by  increasing  the  activation  energy  to  Ea  =  69 A5RT0  for  the 
same  A  used  in  ^1.  The  results  were  nearly  identical  for  all  three 
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Fig.  11.  Reaction  front-propagation  velocities  as  a  function  of  the  position  of  the  leading  edge  of  the  reaction  front  for  x<  300  cm  in  configuration  17.4  calculated  using 
parameter  sets,  (a)  1,  2,  3,  (b)  1,4,  and  (c)  1,  5.  See  Table  4  for  descriptions  of  the  parameter  sets.  Lines  with  symbols  represent  £  =  0.6,  and  those  without  represent  £  =  0.3. 


Table  4 

Parameter  sets  used  in  Figs.  11  and  12  and  their  corresponding  laminar  flame  and  one-dimensional  detonation  properties. 


& 

1 

2 

3 

4 

5 

A  (cm3/g  s) 

1.64  x  1013 

1.2  x  1013 

1.64  x  1013 

4.411  x  1013 

1.64  x  1013 

EJRTo 

67.55 

67.55 

69.45 

67.55 

67.55 

qM/RTo 

39.0 

39.0 

39.0 

34.71 

34.82 

y 

1.197 

1.197 

1.197 

1.197 

1.226 

Si  (cm/s) 

38.02 

32.0 

32.0 

38.02 

38.02 

Tb  (K) 

2210 

2210 

2210 

2000 

2210 

Dq  (m/s) 

1820 

1820 

1820 

1724 

1854 

xd  (cm) 

0.229 

0.331 

0.356 

0.332 

0.0616 
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cases  with  £  =  0.6,  and  only  slight  differences  between  SP3  and  SP\ 
were  found  later  in  the  flame  development  (x  >  250  cm)  for 
£  =  0.3.  These  simulations  suggest  that,  in  the  range  of  5/  consid¬ 
ered,  it  has  negligible  influence  on  the  evolution  of  the  flame  prop¬ 
agation  speed  in  the  channel. 

The  effect  of  the  adiabatic  flame  temperature  on  the  flame 
acceleration  was  tested  by  changing  the  amount  of  heat  released 
per  unit  mass  of  fuel,  while  maintaining  a  constant  5/.  This  was 
done  by  inversely  changing  q  and  A.  Reducing  q  has  the  effect  of 
lowering  both  Sj  and  Tb,  but  a  change  in  A  affects  only  5/.  Thus, 
the  procedure  is  to  set  q  =  (Tb-  T0)/cp  for  the  desired  Tb,  and  then 
to  choose  A  such  that  the  laminar  flame  speed  calculated  using 
the  one-dimensional  flame  model  discussed  in  the  previous  section 
is  equal  to  the  desired  Sb  Parameter  set  SPA  yields  cooler  flames 
(Tb  =  2000  K)  for  the  same  5/ =  38.02  cm/s  as  SP\.  Fig.  lib  shows 
the  flame  acceleration  obtained  using  these  two  parameter  sets 
for  £  =  0.3  and  0.6.  Again,  the  results  obtained  suggest  this  temper¬ 
ature  difference  has  no  noticeable  effect  on  flame  acceleration. 

The  effect  of  varying  y  (and,  hence,  the  sound  speed)  on  the 
flame  acceleration  is  shown  in  Fig.  11c.  Two  sets  of  parameters  that 
give  equal  5/  and  Tb  but  different  ratios  of  specific  heat  are  consid¬ 
ered.  In  the  first  (SPX),  we  use  y  =  1.197.  In  the  second  set,  we  use 
7  =  1.226  and  decrease  q  to  34.82RT0/M  to  maintain  Tb  =  2210  K. 
The  values  of  Ea  and  A  are  unchanged  compared  to  SP\,  so  5/  is  also 
the  same  38.02  cm/s.  Only  slight  differences  between  the  two  cases 
are  noticeable  for  both  f  =  0.3  and  0.6. 

The  second  important  issue  is  how  the  choice  of  model  param¬ 
eters  affects  the  onset  of  the  DDT.  We  note  here  that  it  is  difficult  to 


treat  this  issue  in  a  truly  quantitative  manner  because  DDT  is  a  sto¬ 
chastic  process  that  depends  on  the  formation  of  relatively  small 
hot  spots,  and  these  result  from  combinations  of  shock  reflections 
or  turbulent  fluctuations.  In  some  situations,  even  seemingly 
imperceptible  changes  in  any  physical  or  numerical  parameters 
or  background  conditions  can  lead  to  significant  random  variations 
of  distances  or  times  to  detonation  initiation  [16].  Even  with  these 
caveats,  we  should  still  be  able  to  get  a  qualitative  idea  of  the  like¬ 
lihood  of  DDT  for  a  given  parameter  set.  Fig.  12  shows  the  reaction- 
front  velocities  in  the  downstream  section  of  the  channels  for  the 
cases  shown  in  Fig.  11. 

For  the  cases  shown  in  Fig.  12a,  a  reduction  in  A  or  an  increase 
in  Ea  leads  to  an  increase  in  Xd  (and  X)  and,  hence,  smaller  D  /2  and 
L  //L  When  these  ratios  are  small,  detonations  are  less  likely  to  be 
able  to  propagate  through  the  smaller  space  above  an  obstacle. 
When  a  detonation  fails,  the  reaction  zone  decouples  from  the 
leading  shock  wave  and  propagates  at  a  speed  much  smaller  than 
Dc.  The  length  of  time  the  reaction  front  spends  propagating  as  a 
detonation  is  thus  reduced,  lowering  the  average  propagation 
speed  of  the  reaction  wave.  This  may  account  for  the  slightly  lower 
average  reaction-front  velocities  in  the  quasi-detonation  regime 
(dotted  and  dashed  lines)  shown  in  Fig.  12a  for  f  =  0.3.  A  reduction 
in  A  or  an  increase  in  Ea  also  increases  induction  times  behind 
shocks,  and  this  can  delay  the  detonation  initiation.  Fig.  12a  shows, 
however,  that  the  onset  of  DDT  occurs  later  for  smaller  A,  but  soon¬ 
er  for  larger  Ea.  A  more  comprehensive  parametric  study  would  be 
required  to  determine  conclusively  whether  this  result  is  attribut¬ 
able  to  the  systematic  variation  in  parameters.  For  g  =  0.6,  the  reac- 
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Fig.  12.  Reaction  front-propagation  velocities  as  a  function  of  the  position  of  the  leading  edge  of  the  reaction  front  for  x>300cm  in  configuration  17.4  calculated  using 
parameter  sets,  (a)  1, 2,  3,  (b)  1,4,  and  (c)  1,  5.  See  Table  4  for  descriptions  of  the  parameter  sets.  Lines  with  symbols  represent  £  =  0.6,  and  those  without  represent  £  =  0.3. 
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tion  fronts  propagate  at  the  choking  velocity,  and  all  three  give 
essentially  the  same  propagation  speed. 

Fig.  12b  shows  the  effect  of  simultaneously  decreasing  q  and 
increasing  A.  For  laminar  flames,  this  procedure  lowers  the  adia¬ 
batic  flame  temperature  while  maintaining  constant  5/.  For  detona¬ 
tions,  varying  the  parameters  in  this  way  has  the  effect  of 
decreasing  DCj  and  increasing  xd  and  2.  The  change  in  2  caused  by 
the  change  in  model  parameters  is  not  sufficient  to  drop  D  /2  or 
L  /2  below  their  critical  values  for  the  £  =  0.3  case,  however.  The 
small  effect  on  computed  distances  to  DDT  is  also  consistent  with 
the  fact  that  induction  delays,  which  are  inversely  proportional  to 
qA ,  are  not  significantly  affected.  The  quasi-detonation  propagation 
velocity  for  the  (^4,  £  =  0.3)  case  is  also  less  than  the  result  shown 
for  in  Fig.  12a  for  which  the  values  of  xd  are  similar.  The  differ¬ 
ence  in  the  propagation  velocity  is  likely  caused  by  the  lower  DCj 
attained  using  0>4.  Once  again  the  results  are  similar  for  the 
£  =  0.6  cases,  where  the  reaction  fronts  propagate  at  the  choking 
velocity.  The  reaction  wave  generated  using  ^4  propagates  slightly 
slower  than  that  generated  using  ^1  because  the  temperature 
(and,  hence,  sound  speed)  in  the  combustion  products  is  lower. 

The  choice  of  y  (more  precisely,  y  -  1)  plays  a  large  role  in  det¬ 
onation  initiation  and  propagation.  Using  ^5,  we  investigate  how 
increasing  y  -  1  by  15%  (so  that  y  =  1.226)  affects  the  transition 
to  detonation.  As  discussed  earlier,  in  ^5  a  smaller  value  of  q  is 
used  so  that  the  computed  5/  and  Th  are  the  same  as  those  com¬ 
puted  for  ^1  and  recorded  in  Table  1.  The  average  reaction-front 
velocities  computed  using  ^5  and  0>\  are  shown  in  Fig.  12c.  For 
the  larger-y  simulations,  shock  dynamics  and  post-shock  tempera¬ 
tures  changed  in  the  way  that  reduced  induction  delays  behind 
shocks.  This  is  probably  the  main  reason  why  DDT  occured  sooner 
for  the  £  =  0.3  case,  and  now  even  appeared  for  the  {  =  0.6  case.  The 
computed  value  of  DCj  was  only  about  1.8%  larger  than  that  of  stoi¬ 
chiometric  methane-air  mixtures  (cf.  Table  1 ),  but  xd  was  3.7  times 
smaller  than  that  computed  for  y  =  1.197.  For  £  =  0.3,  the  detona¬ 
tion  propagates  as  a  quasi-detonation  in  a  manner  similar  to  that 
observed  using  set  £P\.  For  f  =  0.6  and  large  y,  we  see  several  in¬ 
stances  of  detonation  initiation  and  failure,  which  is  not  surprising 
since  the  new  D  /2  =  1.5  or  L  /2  =  5.6  are  much  closer  to  their  critical 
values  discussed  in  the  preceding  section.  The  transitions  are  infre¬ 
quent,  and  the  average  propagation  velocity  is  smaller  than  DQ  and 
the  quasi-detonation  propagation  velocities  observed  for  the 
€  =  0.3  cases. 

The  results  shown  in  Figs.  11,  12  indicate  that  varying  the 
model  parameters  to  result  in  relatively  small  (10-15%)  changes 
in  individual  laminar  flame  properties  has  little  impact  on  the  ob¬ 
served  flame  acceleration.  In  general,  the  effects  on  DDT  are  also 
small.  It  is  not  clear  whether  the  differences  in  the  first  occur¬ 
rence  of  DDT  among  the  several  parameter  sets  shown  in 
Fig.  12a  are  due  to  physical  differences  in  the  modeled  systems 
or  to  chance  fluctuations  in  the  thermodynamic  conditions  within 
the  hot  spots  that  initiate  detonations.  The  changes  in  model 
parameters  for  and  both  have  the  effect  of  increasing 
induction  times  behind  shocks  (an  effect  that  could  delay  DDT), 
but  the  resulting  first  occurrences  of  DDT  appeared  later  and 
sooner  than  the  baseline  case,  respectively.  Small  changes  in  the 
thickness  of  the  detonation  wave  (and  detonation  cell  size)  result 
in  correspondingly  small  changes  in  the  average  propagation 
velocity  and  first  occurrences  of  DDT  unless  the  system  is  near 
a  critical  value  for  detonation  propagation.  In  such  a  system,  even 
a  small  increase  in  detonation  cell  size  could  impede  DDT.  In  sys¬ 
tems  sufficiently  larger  or  smaller  than  the  critical  size  for  DDT 
for  a  given  fuel  mixture,  the  interaction  of  the  induced  flow  and 
shocks  with  the  wrinkled  flame  surface  seems  to  have  a  much 
larger  influence  on  the  large-scale  behavior  of  the  reaction  front 
than  the  details  of  the  model  chemistry  in  these  types  of  ob¬ 
structed  channels. 
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There  are  other  physical  properties  of  the  gas  mixture  that  have 
not  been  discussed  here  that  could  potentially  impact  the  first 
occurrence  of  DDT.  For  instance,  the  ratio  of  the  acoustic  time  in 
the  Zeldovich  reactivity  gradient  to  the  reaction  time  [49]  has  been 
shown  to  be  relevant  to  the  development  of  a  detonation  [50].  It  is 
not  possible  to  vary  the  reaction  time  and  the  induction  time  inde¬ 
pendently  using  a  one-step  reaction  model.  This  could  be  done 
with  a  multiple-step  reaction  model,  but  at  a  higher  computational 
cost.  Evaluating  whether  solutions  obtained  using  a  multiple-step 
reaction  model  would  produce  more  accurate  results  and,  if  so, 
whether  the  improvement  justifies  the  additional  computational 
costs  is  left  for  future  investigations. 

4.4.  Flame  resolution  tests 

The  issue  of  how  much  grid  resolution  is  necessary  at  the  reac¬ 
tion  front  becomes  extremely  important  when  attempting  to  sim¬ 
ulate  flame  acceleration  and  DDT  for  very  large-scale  explosions. 
The  size  of  the  system  that  can  be  simulated  will  be  limited  by 
the  number  of  refined  cells  required  to  obtain  a  sufficiently  accu¬ 
rate  description  of  the  reacting  flow.  Thus,  it  is  useful  to  know 
the  largest  cell  size  that  can  be  used  to  meet  this  requirement. 
We  tested  the  effect  of  the  reaction  front  resolution  on  flame  accel¬ 
eration  in  configuration  7.6  with  £  =  1/3  using  the  model  parame¬ 
ters  shown  in  Table  1.  In  all  cases,  the  flame  is  initiated  from  a 
spark  of  size  rf=  0.25  cm.  Three  values  of  dxmin  were  tested, 
0.018124  cm  (corresponding  to  ~3  computational  cells  per  laminar 
flame  thickness),  0.009062  cm  (~5  cells/flame  thickness),  and 
0.004531  cm  (~10  cells/flame  thickness).  Fig.  13a  shows  the  prop- 


(a) 


(b) 


Fig.  13.  (a)  Flame  velocity  and  (b)  flame-surface  length  as  a  function  of  position  of 
leading  edge  of  a  flame  for  configuration  7.6,  £  =  1/3  using  three  successively  more- 
refined  grid  resolutions. 
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agation  velocity  of  the  leading  edge  of  the  flame,  and  Fig.  13b 
shows  the  development  of  the  flame  surface,  both  as  a  function 
of  the  position  of  the  leading  edge  of  the  flame.  As  the  minimum 
grid  size  is  decreased,  and  the  flame  front  becomes  more  resolved, 
the  peak  amount  of  flame  surface  increases.  The  computed  flame 
velocity  during  the  later  stages  of  the  flame  acceleration  process 
increases  as  well  when  the  number  of  grid  cells  in  the  reaction 
zone  is  increased  from  3  to  5.  A  subsequent  increase  from  5  cells 
to  10  cells  has  little  impact  on  the  bulk  acceleration  of  the  flame. 


5.  Summary  and  conclusions 

In  this  work,  we  have  calculated  DDT  in  stoichiometric  mixtures 
of  methane  and  air  in  channels  with  obstacles.  The  systems  simu¬ 
lated  are  too  large  to  compute  with  a  detailed,  multistep,  multispe¬ 
cies  chemical  reaction  mechanism,  and  there  is  no  simple, 
inexpensive  reduced  mechanism  that  can  be  used  for  both  flames 
and  detonations.  Thus,  we  developed  and  calibrated  a  single-step 
reaction-diffusion  model  that  correctly  reproduces  the  length  and 
time  scales  of  both  laminar  flames  and  detonations,  yet  is  compu¬ 
tationally  efficient  enough  to  use  in  large-scale  computations.  The 
computed  velocities  of  the  reaction  fronts  as  the  flames  accelerated 
and  transitioned  to  detonations  showed  reasonable  qualitative  and 
quantitative  agreement  to  those  measured  in  experiments  of  sim¬ 
ilar  sizes  and  geometries.  The  comparison  is  good  enough  to  give 
us  confidence  in  using  the  model. 

Three  stages  of  flame  acceleration  were  observed:  flame 
stretching  and  folding,  flame-front  wrinkling  caused  by  turbulent 
eddies  and  fluid-dynamic  instabilities,  and  flame-surface  creation 
by  shock-flame  interactions.  Higher  levels  of  flame-surface  wrin¬ 
kling  cause  the  flame  to  accelerate  more  rapidly,  but  the  final  chok¬ 
ing  velocity  is  the  same  for  all  geometric  configurations  and 
depends  only  on  properties  of  the  gas  mixture.  Transitions  to  det¬ 
onation  are  observed  in  systems  where  Dorofeev’s  criterium 
r  =  (S  +  D)/2/(l  -  D*/D)  was  large  compared  to  the  detonation  cell 
size.  Detonations  appear  when  shock  reflections  from  channel 
walls  and  obstacles  locally  raise  the  temperature  in  the  unburned 
gases  to  the  ignition  point.  This  initiates  a  reaction  that  can  ignite 
a  detonation  by  Zeldovich’s  gradient  mechanism.  The  survival  of 
the  detonation  then  depends  on  the  local  thermodynamic  condi¬ 
tions  and  the  size  of  the  orifice  through  which  it  will  propagate. 
Frequent  detonation  failure  and  reignition  lowers  the  average  det¬ 
onation  velocities  to  below  DCj,  which  is  consistent  with  experi¬ 
mental  evidence. 

We  showed  that  the  large-scale  flame  dynamics  and  DDT  are 
generally  insensitive  to  small  changes  in  model  parameters.  Only 
relatively  minor  differences  in  the  acceleration  of  the  flames 
were  observed  for  10-15%  variations  in  the  laminar  flame  speed, 
adiabatic  flame  temperature,  and  specific-heat  ratio.  The  dis¬ 
tances  to  DDT  showed  a  somewhat  wider  variation  for  different 
parameter  sets.  The  most  likely  explanation  for  this  is  the 
dependence  of  ignition  delay  times  behind  shocks  on  the  model 
parameters.  Reflected  shocks  and  Mach  stems  are  more  likely  to 
be  able  to  ignite  a  detonation  when  the  ignition  delays  are  short, 
but  whether  or  not  this  occurs  at  a  particular  location  is  highly 
dependent  on  the  system  geometry  and  stochastic  variations  of 
local  thermodynamic  conditions.  This  makes  the  exact  prediction 
of  the  location  of  DDT  extremely  difficult.  Definitively  separating 
the  effect  of  the  model  parameters  on  DDT  from  the  effects  of 
stochasticity  would  require  a  more  comprehensive  study  and  is 
left  for  future  work. 

The  parameter  that  has  the  most  pronounced  effect  on  DDT  is 
the  specific  heat  ratio.  The  value  of  y  determines  post-shock  tem¬ 
peratures  and  pressures,  which  strongly  influence  the  detonation 


of  a  ID  detonation  by  a  factor  of  3.  For  the  systems  considered,  this 
change  affected  the  distance  to  DDT  for  £  =  0.3  and  the  ability  of  a 
detonation  to  form  when  £  =  0.6.  The  strong  dependence  of  detona¬ 
tion  structure  on  y  makes  the  choice  of  this  parameter  important 
for  the  simulations.  A  future  area  of  investigation  would  be  to 
see  how  using  a  temperature-dependent  y  would  affect  shock- 
flame  interactions  and  DDT. 

A  series  of  grid  resolution  tests  showed  that  increasing  the  res¬ 
olution  from  5  to  10  cells  per  laminar  flame  thickness  has  practi¬ 
cally  no  effect  on  the  flame  acceleration.  Small  differences  in  the 
flame  speed  found  for  under-resolved  reaction  zones  (3  cells  per 
laminar  flame  thickness)  have  only  a  minor  impact  on  the  large- 
scale  flame  development. 

In  spite  of  geometrical  differences  between  the  model  configu¬ 
rations  used  in  the  simulations  and  the  experimental  systems  to 
which  they  are  compared,  we  observed  reasonable  agreement  in 
flame  development  and  DDT.  Such  quantitative  agreement  be¬ 
tween  2D  simulations  and  3D  experiments  may  seem  strange  since 
the  difference  in  the  equilibrium  energy  cascade  for  two-dimen¬ 
sional  and  three-dimensional  turbulence  is  substantial.  It  has  been 
suggested,  however,  that  in  flows  where  repeated  shock-flame 
interactions  drive  the  fluid  instabilities,  the  turbulence  is  nonequi¬ 
librium  and  that  energy  is  transferred  directly  into  a  broad  range  of 
scales  simultaneously  (see,  e.g.,  [15]).  Under  this  paradigm,  differ¬ 
ent  scales  of  turbulent  motion  are  populated  more  efficiently  than 
through  the  standard  Kolmogorov  energy  cascade  model.  In  our 
simulations,  shock-flame  interactions  become  important  as  the 
velocity  of  the  leading  edge  of  the  flame  approaches  the  speed  of 
sound  in  the  burned  material  and  play  a  large  role  in  the  final  stage 
of  the  flame  acceleration  process. 

DDT  in  our  simulations  occurs  due  to  the  appearance  of  hot 
spots  behind  a  shock  wave  or  Mach  stem.  Whether  or  not  these 
hot  spots  can  trigger  a  detonation  depends  on  the  strength  of  the 
shock.  We  found  that  the  critical  system  sizes  for  DDT  in  our  calcu¬ 
lations  were  smaller  than  those  in  the  experiments.  The  level  of 
uncertainty  inherent  in  the  determination  of  average  detonation 
cell  sizes  in  methane-air  mixtures,  however,  precludes  us  from 
making  a  quantitative  assessment  of  DDT  criteria. 
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Appendix  A.  ID  Steady-state  flame 

We  assume  that  the  Lewis  number  of  the  fuel  mixture  is  unity, 
which  gives  Y  as  a  linear  function  of  temperature, 

Y  =  l-(T-T0)Cp/q,  (A.l) 

so  that  the  the  reaction  rate  becomes  a  function  of  temperature 
only,  Q  =  A(1  -  (T  -  T0)Cp/q)exp(-Ea/RT). 

We  change  the  boundary  value  problem  (Eqs.  ( 1 1 ),  ( 1 2))  into  an 
initial  value  problem  by  solving  the  inert  heat  transport  equation 


structure.  A  15%  change  in  y  -  1  was  found  to  decrease  the  width  (derived  by  neglecting  Q  in  Eq.  (11)), 
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r)F 

ffT  =  Pucr’  (A-2) 

for  T  close  to  T0  to  obtain  the  temperature  gradient  in  the  vicinity  of 
T0  for  an  initial  guess  of  S/.  Eq.  (A.2)  is  integrated  from  T0  to  T  using 
a  second-order  Runge-Kutta  method  with  variable  step  AT  that 
exponentially  increases  with  T.  The  distance  x  is  integrated  simulta¬ 
neously  using  Eq.  (12).  The  choice  of  T  has  little  effect  on  the  final 
solution  as  long  as  qpT2(T)  is  small  compared  to  DCpdT/dx. 

The  computed  x  and  dT/dx  at  T  are  used  as  initial  conditions  to 
Eqs.  (11)  and  (12),  which  are  integrated  using  a  second-order  Run¬ 
ge-Kutta  method  with  constant  step  Ax.  The  temperature  gradient 
at  T=Tb  =  T0  +  q/Cp  depends  on  Sh  and  the  integration  procedure  is 
repeated  with  different  Si  until  the  condition  dT/dx  =  0  when  T  =  Tb. 

Appendix  B.  ID  Steady-state  detonation 


For  a  ID  planar  shock  wave  moving  at  DCj  (Eq.  (16)), 


pznd  _-,M2  y  y  - 1 
p0  gr  +  i  y  + 1’ 

(B.l) 

Pznd  _  Mcj  (7  +  1) 

Po  M2CJ(y- 1)  +  2’ 

(B.2) 

&znd  —  —  0.5 (Pznd  +  Po )  ( 1  / p0  ~  1/ Pznd )  ■ 

(B.3) 

where  MC;  =  DC;/c0  and  e0  =  P0/(p0(y  -  1)).  These  ZND  parameters 
( PZnd ,  Pznd ,  and  eZND)  are  used  as  initial  conditions  for  the  integra¬ 
tion  of  Eqs.  (13)— (15).  The  equations  are  integrated  using  a  sec¬ 
ond-order  Runge-Kutta  method  from  the  initial  conditions  to  the 
position  x*  where  U  =  DQp0lp  =  c.  The  half-reaction  thickness,  xd  is 
then  defined  as  the  distance  between  the  leading  shock  wave  and 
the  point  where  half  of  the  fuel  has  been  consumed  in  the  flame 
zone,  i.e.,  where  Y  =  0.5. 
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